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ABSTRACT 

It has been known for some time that rotating bars in galaxies slow due to dynamical friction against 
the halo. However, recent attempts to use this process to place constraints on the dark matter density 
in galaxies and possibly also to drive dark matter out of the center have been challenged. This paper 
uses simplified numerical experiments to clarify several aspects of the friction mechanism. I explicitly 
demonstrate the Chandrasekhar scaling of the friction force with bar mass, halo density, and halo velocity 
dispersion. I present direct evidence that exchanges between the bar and halo orbits at major resonances 
are responsible for friction and study both individual orbits and the net changes at these resonances. I 
also show that friction alters the phase space density of particles in the vicinity of a major resonance, 
which is the reason the magnitude of the friction force depends on the prior evolution. I demonstrate that 
bar slow down can be captured correctly in simulations having modest spatial resolution and practicable 
numbers of particles. Subsequent papers in this series delineate the dark matter density that can be 
tolerated in halos of different density profiles. 

Subject headings: galaxies: formation — galaxies: kinematics and dynamics — galaxies: halos — dark 
matter 



1. INTRODUCTION 

Friction between a rotating bar and a massive halo was 
first reported many years ago (Sellwood 1980) and has 
subsequently been worked on sporadically (Tremaine & 
Weinberg 1984; Weinberg 1985; Little & Carlberg 1991; 
Hernquist & Weinberg 1992; Athanassoula 1996). It has, 
however, received a lot of attention in recent years as a 
potential probe of, and structuring mechanism for, dark 
matter halos. Debattista & Sellwood (1998, 2000) argued 
that the fact that bars appear not to have been slowed 
places an upper bound on the density of the dark matter 
halo in barred disk galaxies. Weinberg & Katz (2002) 
argue that the transfer of angular momentum from the 
bar to the halo could reduce the central density of the 
dark matter halo by a substantial factor. 

Both of these claims have subsequently been disputed; 
Valenzuela & Klypin (2003) claim a counter-example of a 
bar that does not experience much friction in a "cosmo- 
logically- motivated" halo and HoUey-Bockelmann & Wein- 
berg (2005) find that cleverly- constructed uniform-density 
halos can also avoid friction. The argument by Athanas- 
soula (2003) that weak bars experience little friction is 
clearly correct, but irrelevant for strong bars. Only mod- 
est reductions in halo density have been achieved so far 
(HoUey-Bockelmann, Weinberg & Katz 2003) and Sell- 
wood (2003) found that contraction of the disk mass distri- 
bution as it lost angular momentum to the halo dragged 
halo mass in with it, causing the density to rise, rather 
than to decrease. 

Valenzuela & Klypin (2003) and Holley-Bockelmann, et 
al. blame discrepant conclusions on inadequacies respec- 
tively of the codes and of the number of particles employed, 
but it is also likely that a good part of the differences re- 
ported by these authors arise because the physical models 
differ. There have been a few tests with different codes 
from the same initial conditions and the results compare 
quite well {e.g. O'Neill & Dubinski 2003; Sellwood 2003); 



other comparisons are reported here and in Papers II & 
HI (Sellwood & Debattista 2005a,b). 

The present paper attempts to clarify the physics of dy- 
namical friction between a bar and a halo and to show 
that it can be correctly reproduced in simulations of a 
size that is readily accessible with current computational 
resources. The counter-example reported by Valenzuela 
& Klypin (2003) is addressed in Sellwood & Debattista 
(2005c) and Paper II. The constraint on halo density that 
can be deduced from the existence of strong, fast bars and 
criticisms by Athanassoula (2003) are addressed in Paper 
HI. 

2. THEORETICAL BACKGROUND 

Dynamical friction (Chandrasekhar 1943) is the re- 
tarding force experienced by a massive perturber moving 
through a background of low-mass particles. It arises, even 
in a perfectly coUisionless system, from the vector sum of 
the impulses the perturber receives from the particles as 
they are deflected by its gravitational field. Equivalently, 
friction can be viewed as the gravitational attraction on 
the perturber of the density excess, or wake, that develops 
behind it as it moves, as was nicely illustrated by Mulder 
(1983). 

Chandrasekhar's formula (eq. 7-17 of Binney & Tremaine 
1987; hereafter BT) for the acceleration qm of a per- 
turber of mass M moving at speed vm through a uniform 
background, density p, of non-interacting particles having 
an isotropic velocity distribution with a 1-D rms velocity 
spread a, may be written as 



= 47rlnAG2^yr^ 



(1) 



Here, A {— bmax/^min) is the argument of the usual 
Coulomb logarithm, and the dimensionless function V is 
drawn in Figure ^ for a Gaussian distribution of veloc- 
ities; other velocity distributions would yield a different 
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Fig. 1. — The dimensionless acceleration function V defined in 
equation Q for the case of a Gaussian distribution of velocities 
among the background particles. 

functional form. 

The simplifying assumptions in its derivation strictly in- 
validate application of equation to the physically more 
interesting problem of friction in a non-uniform medium in 
which the background particles are confined by a potential 
well and interact with the perturber repeatedly. Repeated 
encounters between the perturber and the background par- 
ticles require a more sophisticated treatment. 

2.1. Orbits and resonances 

The motion of a particle pursuing a regular orbit in 
a smooth ellipsoidal potential is triply periodic, and is 
most conveniently described by action-angle variables (see 
BT).^ The three angular momentum-like actions, J, de- 
scribe the amount of "round and round" , "in and out" and 
"up and down" motion associated with the orbit, while 
the angles, w, specify the phases of these three separate 
oscillations that underlie the motion. One of the many 
advantages of these variables is that the angles are defined 
such that they increase at a uniform rate w = 0( J) for 
each orbit. 

In this paper, as in most other discussions, we will be 
concerned with mild perturbations to a spherical potential, 
where the unperturbed motion of any particle is confined 
to a plane. In this case, there are just two non-zero actions: 
the azimuthal action = L, the total specific angular 
momentum, and the radial action, J^. The motion is then 
simply doubly periodic at angular rates and fJ^- 

The angular frequency, Q^, of a particle's mean motion 
about the center is, in general, incommensurable with its 
angular frequency of radial motion, fir- However, all or- 
bits can appear to close when viewed from many different 
rotating frames. An observer rotating at the rate 

n' + kQr/m, (2) 

would see the orbit close after m radial oscillations and k 

^Irregular orbits that do not preserve three integrals of motion 
may be important in real galaxies, but most discussion of dynamical 
friction has focused on simple (usually spherical) potentials in which 
all unperturbed orbits are regular. 



turns about the center; Kalnajs (1977) draws an orbit in 
several of the most important frames. 

The angular momentum vectors of orbits in a spherical 
system are distributed over all angles. We define some 
direction as the ^-axis and use 6 for the angle between 
that axis and the plane of each unperturbed orbit, with 9 
having the same sign as the z-component of the angular 
momentum Lz = LsinO. The frequencies of the motion 
projected into the equatorial plane, 9 = 7r/2, are indepen- 
dent of 9, save only that the sign of fl^ is that of L^. 

We are interested in an ellipsoidal bar-like perturbation 
that rotates about its minor axis, which we take to be the 
z-axis. If the bar rotates at an angular speed Qbi there 
are many possible resonances between the orbits of the 
particles and the perturbing potential. Resonances occur 
where mf2{, = nO^ -I- k^r, where (m, k,n) = m are three 
integers; m is even for a bar, n = 0, ±to, and — 1 < A: < oo 
(Tremaine & Weinberg 1984). Three of the most impor- 
tant resonances are familiar from disk dynamics: the coro- 
tation resonance (CR) where n = ni and fc = 0, the inner 
Lindblad resonance (ILR) where n = m = 2 and fc = — 1, 
and the outer Lindblad resonance (OLR) where n = m = 2 
and k = 1. Negative values of n lead to resonances with 
retrograde particles, which seem to be of little dynamical 
importance. Weinberg & Katz (2005; hereafter WK05) 
stress the importance a fourth: the direct radial resonance 
(DRR) for which n — and k = 1. This last resonance is 
absent for orbits in the bar plane and is most important 
for near polar orbits. 

2.2. Friction in spheroidal systems 

Tremaine & Weinberg (1984, hereafter TW84) devel- 
oped the basic theory of dynamical friction by a generic 
perturber in spherical systems. Following the precepts of 
Lynden-Bell & Kalnajs (1972), they derived an expression 
for the "LBK torque" in a spherical system, which they 
showed arises purely at resonances. Weinberg (1985) gave 
a specific evaluation of the LBK torque for the case of a 
bar. The role of resonances has been re-emphasized in 
recent work (Weinberg & Katz 2002; Athanassoula 2003; 
HoUey-Bockelmann, et al. 2003; WK05). 

As for the classical Chandrasekhar problem, friction can 
be viewed in two equivalent ways: it can be seen either as 
arising from the gravitational coupling between the bar 
and a misaligned density response in the halo, or it can be 
interpreted as the net effect of exchanges at resonances be- 
tween the particles and the perturber. The misalignment 
can be understood in terms of interactions at resonances, 
where halo particles gain or lose angular momentum. 

The daunting expression for the LBK torque derived 
by TW84 (their eq. 65) does not appear to resemble 
eq. iQJ above. However, TW84 also reformulate the orig- 
inal Chandrasekhar problem in a manner that highlights 
the similarities with the LBK torque in a spherical system, 
but which contains no resonant terms because the system 
is infinite. A nai've guess, therefore, is that the LBK torque 
causes the angular acceleration of the perturber to scale 
as 

^ MnPief^)^ (3) 

where Alp is the mass of the perturber, ps and are the 
characteristic density and velocity dispersion of the spher- 
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ical system, and the dimensionless function Q contains all 
the complicated dependence on the details of the distribu- 
tion function, potential well, resonances, bar shape, etc. A 
characteristic length scale a is included in order to make 
the argument of Q dimensionless. The fmictions Q{x) and 
V{x) share the general properties that they are negative, 
at least for isotropic distribution functions (TW84), they 
must ^ as a; ^ c», and that they can reasonably be 
expected to be oc a; as a; ^ 0. 

Weinberg (2004) points out that the LBK torque for- 
mula derived by TW84 is not the full story, because their 
expression is evaluated for a perturbation of fixed ampli- 
tude and pattern speed. Weinberg finds that the torque 
depends strongly on the previous time-dependence of the 
perturbation, which he needed to take into account in or- 
der to reconcile his perturbation theory with results from 
his simulations. This improvement to the theory is a major 
step forward, as we show here. 

2.3. Motivation 

The numerical experiments of Lin & Tremaine (1983) 
showed that a satellite orbiting a spherical system com- 
posed of much lighter particles confined in a potential well 
experiences a frictional force that scales very much as pre- 
dicted by eq. or eq. Q, despite the complications 
caused by resonances. This, and other evidence, led BT 
to conclude that "Chandrasekhar's formula often provides 
a remarkably accurate description of the drag experienced 
by a body orbiting in a stellar system" . 

We might hope that the Chandrasekhar scaling would 
also hold for bars, but there is no detailed check in the lit- 
erature. Weinberg (1985) reports a few crude simulations 
of the Lin & Tremaine type in support of his theoretical 
calculations and he remarks (Weinberg 2004) that he has 
recently made more such tests. 

In §Siynil present a much more extensive study using 
the Lin & Tremaine technique for this slightly different 
physical problem, and show that the expected scaling of 
eq. © holds quite well. reports a more detailed study 
of resonant exchanges both of indiviual orbits and of the 
ensemble of particles. fQconfirms Weinberg's (2004) find- 
ing that the time dependence of the perturbation is of 
major importance, and offers an explanation. I present 
convergence tests without self-gravity in fjSl and I include 
the effects of self-gravity between the halo particles in 
in order to determine empirically the numbers of particles 
needed to reproduce the correct frictional force in fully 
self-consistent simulations. 

3. MODELS AND METHOD 

3.1. Halo models 

I create an iV-body realization of a spherical mass dis- 
tribution, which for brevity I describe as a halo, although 
it could be any spherical system of coUisionless particles. 
In 21 the particles move in the smooth analytic gravi- 
tational potential of the adopted halo and I neglect any 
interaction forces between the particles. Since I draw the 
particles from a distribution function (DF) that generates 
the adopted halo density, the particle distribution is in 
equilibrium and does not evolve in the absence of an ex- 
ternal perturbation. The isotropic halos start with no net 
angular momentum. 



I use three quite different halo models. The first is a 
Hernquist (1990) model, which has the density profile 
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(4) 



with total mass Mh- The density profile declines as 
for r ^ th and as r~'* for r ^ th . It should be noted that 
this model differs only slightly from the NFW profile (see 
Appendix), which appeared to be a reasonable fit to early 
simulations (Navarro, Frenk & White 1997) of the collapse 
of dark matter halos. Hernquist gives the expression for 
the isotropic DF that generates the halo of eq. I^J. I do not 
employ the infinite model, but remove from the DF any 
particle with sufficient energy to reach r > 20rH, so that 
the active mass in particles is ~ 0.86Mh while they move 
in the analytic potential of the untruncated halo. The 
active density profile is very little affected for r < ISrjj, 
while the bars I employ are typically much smaller, with 
semi-major axes ~ rpj. 

I have also employed the well-known Plummer sphere 
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(5) 



which has a uniform density core. The isotropic DF that 
generates this density profile is a polytrope of index n = 5 
(BT, equation 4-104). Again I eliminate any particle with 
sufficient energy to reach r > 20rp, which is less than 1% 
of the mass. 

As a third halo model, I have adopted the singular 
isothermal sphere (SIS), which formally has the scale- free 
density profile 
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An isotropic DF that generates this density has a 1-D ve- 
locity dispersion a ~ Vq/ \f2. The particles move in the ex- 
act logarithmic potential of the untruncated sphere, but I 
generate active particles from this DF with a limited range 
of energies. The upper bound is set so that no particle 
has enough energy to pass an outer radius Ti^axi while the 
lower bound eliminates any particle that would be bound 
inside some small radius Tmin- I choose Tmax = 20a and 
''min = O.Ola, where a is my adopted bar semi-major axis. 
For this model only, I set the mass of each halo particle 
proportional to L^/^, where L is its total angular momen- 
tum, in order to obtain a disproportionately higher density 
of particles in the inner parts of the halo. 

3.2. Bar 

The bar model is a homogeneous ellipsoid, which has 
the density distribution 
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where M{, is the mass of the ellipsoid, a, b c are its 
semi-axes, with a > b > c, and 
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I do not use the full field of this bar but employ only non- 
axisymmetric parts of the field, excluding the monopole 
term in order not to disturb the radial profile of the halo 
as I introduce the bar. This approach is similar to that first 
adopted by Hernquist & Weinberg (1992), who employed 
only an approximation to the quadrupole term of the bar 
field. Here, I determine the precise field of the bar using a 
multipole expansion {e.g. BT, §2.4), but use only the non- 
axisymmetric quadrupole (1 — 2, m = 2) term, and in 
some cases higher terms, to accelerate the halo particles. 
(The odd-Z and odd-m terms all vanish because the bar is 
respectively symmetric about the mid-plane and has 2-fold 
rotational symmetry, while the (/, 0) terms do not rotate. 
Because the bar lies in the equatorial plane, terms with 
/ — m will be much larger than those with I > m > 0.) 
The internal and external contributions to the bar field 
over the radial range of the bar are tabulated only once 
and stored; it is straightforward to rotate the tabulated 
values through any desired angle at each step. 

The bar rotates at angular frequency fif, about its short- 
est axis, and I introduce the bar smoothly by increasing 
Mb as a cubic function of time from zero at t = to its 
final value a,t t = tg. I consider the model bar to be rigid, 
and to spin down due to loss of angular momentum ac- 
cording to its moment of inertia about the shortest axis, 
which is 

I=^{a' + b'). (9) 

A rigid bar is essential for this study, even though it is 
unrealistic in many ways. I must make arbitrary choices 
for the bar mass, length, axis ratio, density profile, and 
initial pattern speed, rather than have all these parameters 
set by the dynamics of a disk. In addition, a rigid bar will 
not have the same moment of inertia as a figure of the same 
density profile composed of active particles, which has a 
pattern speed set by the mean precession rate of the bar 
particle orbits. Furthermore, the density profile of the bar 
should change as the bar loses angular momentum and the 
mean radius of the bar particles decreases, and a bar in an 
active disk may also grow by trapping extra orbits. The 
present study, however, requires rigid bars with controlled 
parameters in order to determine how friction depends on 
the bar properties. 

3.3. Numerical procedure 

The halo particles move in the combined fields of the 
halo and of some non-axisymmetric terms of the rotating 
bar field. The halo field is rigid in the experiments de- 
scribed in 1^ but some self-gravity terms are included in 
^ I sum the z-component of the net torque on the parti- 
cles and use the negative of this sum as the torque acting 
to accelerate the bar at every step, so that the combined 
angular momentum of the halo and bar is conserved. 

My system of units is such that G = Mh = = 1, 
where rx = for the Hernquist halo and r^ = rp for the 
Plummer halo. For the SIS models, on the other hand, 
I choose G = Vq = a = 1, for which my unit of mass is 
qVq /G. Unless otherwise stated, the bar axes are a : b : 
c — 1 : 0.5 : 0.05 and tg = 10 in these units. 

As particles in these models have a wide range of fre- 
quencies, I divide the computation volume into a number 
of spherical zones (typically 5) in which particles move 
with different time steps that differ by factors of 2. The 
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Fig. 2. — (a) The variation of with angular speed of the bar. 
The sohd and dotted curves show the acceleration in response to 
the quadrupole field of a fat (1:0.5:0.05) bar for different initial flf, 
values. The dashed curve shows the same quantity when the per- 
turbing force is the (4,4) component of the bar field only, (b) The 
phase lag in radians between the density response in the inner halo 
and the bar from the run shown by the solid curve in (a). The re- 
sponse is approximately orthogonal to the bar direction when Qi, is 
large, and gradually shifts into alignment as the bar slows. 

integration step changes when particles cross zone bound- 
aries. I have verified that the results presented here do not 
depend on this scheme, or on the radii of the zone bound- 
aries, or on the adopted fundamental time step, neither are 
they affected when the bar quadrupole field is replaced by 
the smooth function adopted by Hernquist & Weinberg 
(1992). 

4. EXPLORATION OF PARAMETER SPACE 
4.1. Scaling with angular speed 

The solid curve in Figure [S^a) shows the function 8 for 
a Hernquist halo and a bar with axis ratios a : b : c = 1 : 
0.5 : 0.1 and a — rjj, when Qh — 1.5 initially. (The bar 
mass Mf, = O.OlMh and the halo is represented by lOM 
particles.) Over most of the range, the drag force on the 
bar varies quite smoothly, peaking when flf, ~ 0.8. An 
initial transient, associated with the turn-on of the bar, is 
evident, and the curve also has a feature near fib ~ 0.2. 

In this case, the initial pattern speed of 17;, = 1.5 is 
unrealistically high, and corotation lies well inside the bar, 
at r ~ 0.274a. The peak deceleration of the bar occurs 
when corotation lies at the still unrealistically small radius 
r ~ 0.58a. The pattern speed must drop to ilb ~ 0.5 to 
place corotation at the end of the bar, by which time the 
drag force is roughly one third of its peak value. 

The dotted curves in Fig. HJa) show, for identically the 
same bar and halo, the acceleration in a number of other 
experiments with lower initial starting speeds. The sur- 
prisingly large differences confirm Weinberg's (2004) find- 
ing that the friction force depends on the past evolution. 
The largest differences arise at angular speeds high enough 
for corotation to lie inside the bar and differences are 
smaller when fib <, 0.5. I discuss this behavior more fully 
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Fig. 3. — As for Fig. I^a) but for a thinner bar with axis ratios 
1:0.2:0.05. Note the different scales between the two figures. 

in 371 

Fig. Olb) shows the lag angle between the principal axes 
of the bar and of the quadrupole response in the halo as 
a function of Qb for the case when flf, = 1.5 initially. I 
estimate the position angle of the halo density response 
from the phase of the (2,2) component of a high radial or- 
der spherical Bessel function {e.g. Arfken 1985) transform 
of the particle distribution. It is clear that the response 
lags the bar by almost a right-angle at high fif,, and the 
lag angle generally decreases as the bar slows, becoming 
nearly aligned with the bar when its rotation speed is very 
low. Comparison with Fig. |2Ia) confirms, as it must, that 
the torque is weak when the response is almost orthogonal 
to, or aligned with, the bar and is greatest as the phase 
lag passes through intermediate angles. (Since the magni- 
tude of the torque also depends on the amplitude of the 
response, the maximum need not be when the response is 
precisely 45° out of phase, although it clearly occurs close 
to this angle.) 

The phase angle of the halo response shown in Fig. |2Ib) 
is not hard to understand. The following argument is spe- 
cific to spherical potentials, but it generalizes to aspherical 
cases. 

As already noted, any unperturbed orbit will close in 
a frame that rotates at the rate Vtm = {nfl^ + kQ.r)/m 
(n ^ 0), but only those orbits for which Vim = also 
close in the frame of the rotating perturbation, and are 
exactly in resonance. 

Most orbits are not resonant, however, and ^ for 
any to. For these general orbits, the perturbation adds to 
its otherwise axisymmetric time-averaged density, a forced 
non-axisymmetric distortion that corotates with the bar. 
As happens for simple harmonic oscillators, the driven re- 
sponse is in phase, or aligned, with the bar when f2(, < f2mj 
and is perpendicular to the bar for higher bar pattern 
speeds. Thus the forced response of an orbit switches from 
perpendicular to alignment as the bar slows across its res- 
onant frequency; the change of phase is gradual because 
the resonance is broadened by the changing pattern speed. 
Since many orbits are present with a wide range of pre- 




FlG. 4. — As for Fig. 13 a), but (a) for the Plummer halo and (b) 
for the singular isothermal sphere. 

cession rates for each resonance, the net response is the 
aggregate of many orbits. 

4.2. Different strength bar 

Fig.Oshows the much stronger frictional deceleration for 
a thinner bar with axis ratios 1:0.2:0.05. The quadrupole 
potential of this bar peaks at about twice the value of 
the fatter bar used in Fig. Ola), but at a smaller radius; 
the two fields could not be matched by scaling, therefore. 
Nevertheless, the curve has a similar shape. 

4.3. Four-fold and higher potential terms 

The dashed curves in Figs.|2Ia) &|31show the frictional 
deceleration for a 20 times more massive bar when only 
the Z = 4, 771 = 4 term of the bar field is used to force the 
halo. Since evolution with a 1% mass bar would be im- 
practicably slow, I used a 20% mass bar in both cases, and 
scaled the curves appropriately. The relative strengths of 
the quadrupole and higher-order components of the per- 
turbing field depend on the bar shape; the (4,4) component 
of the potential has a peak amplitude of ~ 1/6 that of the 
quadrupole for the fatter bar, whereas it is a larger frac- 
tion 1/3) for the thinner bar. Yet the contribution to 
friction from the (4,4) term is small even for the sharper 
bar (Fig.O; higher order terms are still less important. 

Since friction is dominated by only the lowest-order, 
non-axisymmetric term of the potential, it is possible for 
even low-spatial resolution codes, such as that used by 
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Fig. 5. — The bar pattern speed as a function of normalized 
time for a series of experiments with differing mass bars, but all 
other properties held fixed. Bar masses are 0.5%, 1%, & 2% (all 
solid lines) and 5% (dotted line) of the halo mass. The times are 
scaled by the bar masses and shifted horizontally to coincide when 
Cfc = 0.3. 

Debattista & Sellwood (1998, 2000), to reproduce the 
frictional drag quite accurately. Thus the suggestion by 
Valenzuela & Klypin (2003) that their different result was 
due to inadequate spatial resolution in the earlier work is 
unlikely to be correct. The different results obtained in 
these two studies will be discussed further in Paper II. 

4.4. Other halo density profiles 

Figure 2J a) shows the variation of with angular speed 
for the Plummer halo, while Fig. 0{b) is for the singu- 
lar isothermal sphere (SIS). These are directly compara- 
ble with Fig. [2Ia), which is for the Hernquist halo. The 
angular velocity dependence of the frictional acceleration 
is noticeably different in both cases. 

The bar in the Plummer halo has a semi-major axis 
a — Rp and a mass Mf, — Q.02Mh. Friction has two ap- 
proximately equal peaks when fif, ~ 1.4 and Qb — 0.7, but 
the angular speed of this bar must drop to 2^'^/'* ~ 0.6 
before corotation moves outside the bar. 

The behavior for the SIS is different again. Corotation 
is at the end of the bar when fif, = 1, where friction is 
again past its peak. Note that features can arise in this 
curve because the scale-free nature of this model is broken 
in two ways: the perturbation has a definite linear size and 
the energy range of the active particles is restricted. The 
mass of the bar in this case. Mi, = 0.2aV^^/G. 

These curves, and that in Fig. [S^a), indicate the func- 
tional form of Q in eq. © . Since the bar is the same in all 
three cases, the differences stem directly from differences 
in the halo mass profiles, which affect the frequencies and 
particle density at the different resonances. Such differ- 
ences are the analog of having a different velocity distri- 
bution for the background particles in the Chandrasekhar 
problem. 

5. SCALING WITH OTHER PARAMETERS 

In this section, I check directly whether the accelera- 
tion scales with bar mass, halo density and halo velocity 
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Fig. 6. — As for Fig. |21 but for the case when the halo density 
is scaled. The dotted curve is for the case when the halo density is 
doubled. See i|5.2l for details. 

dispersion, as suggested by eq. ©. All experiments re- 
ported in this section employ lOM particles, as were also 
used for those shown in Figs. [3 El but the initial bar 
pattern speed is set at the more conventional value such 
that corotation is at the end of the bar. I illustrate these 
tests for the Hernquist halo; the other halo models scale 
about as well, except for the test with velocity dispersion, 
as discussed in ^5.'6\ 

5.1. Scaling with bar mass 

Figure [S] shows the time evolution of the pattern speed 
for bars of different masses, in the Hernquist halo; the bar 
pattern speed is a smoother function than its acceleration. 
Each line in this Figure is from a run with a different bar 
mass in the range 0.005 < Mt/Mh < 0.05 and fib = 0.5 
initially; the time axis is scaled by Mb, and the curves have 
all been shifted horizontally to coincide at the moment 
at which fib = 0.3. The curves overlay almost perfectly, 
indicating that frictional acceleration scales linearly with 
Mb in this regime, as expected from eq. (j^J. 

5.2. Scaling with halo density 

Equation ^ predicts that the acceleration is propor- 
tional to the density of the background. In the present 
more realistic case, the density of the background is a 
function of radius and it might seem that I would have 
to evaluate the LBK torque for a number of different ha- 
los to obtain a testable prediction. Fortunately, there is 
a much simpler way to test for the density dependence: I 
simply scale the density of the entire halo without chang- 
ing the potential well in which the particles move. This 
ploy is equivalent to treating some fraction of the halo as 
rigid if the density is reduced but, since the density and 
potential of the halo need not be self-consistent in these 
restricted experiments, it is also possible to increase its 
density. I adopt this admittedly artificial strategy here 
simply to test for the expected linear scaling. 

Figure El shows the results with a 2% mass bar, with 
Hernquist halos that have 0.5, 1 and 2 times the density 
given in eq. Q. Again the time axis scale is proportional 
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Fig. 7. — (a) As for Fig.|^ but for the case when the halo velocity 
dispersion and potential are scaled. The dotted curve is for the 
largest velocity dispersi on, (b) Same as (a), but for the singular 
isothermal sphere. See t|5.3l for details. 



to the density, and curves are shifted so that they coincide 
when fif, = 0.3. While not quite as convincing as Fig. \S\ 
the similarity of the curves shows that the expected lin- 
ear scaling holds approximately, in satisfactory agreement 
with the prediction of eq. (PJ. 

5.3. Scaling with halo velocity dispersion 

Finally, I test the scaling with velocity dispersion pre- 
dicted by eq. l|2Jl. Again I employ a trick: by scahng the 
gravitational potential by a factor, a^, and the velocity 
dispersion, cr, by a, we continue to have an equilibrium 
model with an unchanged density profile. (The halo den- 
sity and potential are no longer self-consistent, of course.) 
For this series of tests, however, I must also scale the bar 
pattern speed, since the argument of Q contains the dis- 
persion, a. 

Figure [TJa), which plots fib/a as a function of Mta~^, 
shows that when the velocity dispersion is increased or de- 
creased by a factor \/2, the evolution of the scaled bar an- 
gular speed does vary at approximately half or double the 
rate, respectively, as predicted by eq. ||2J), but the curves 



do not match up as well as in Figs.|Sl&|Sl 

Since scaling the potential well affects all the resonances, 
the particle density at each resonance does not scale in a 
simple manner, and formula Q is too naive. However, 
scaling should be restored in a scale-free model. The same 
scaling test with the SIS supports this expectation; the 
curves in Fig. |7Jb) follow each other more closely than in 
Fig.[7fa), which was for the Hernquist halo. 

5.4. Summary 

The present problem differs from that considered by 
Chandrasekhar in many ways: the perturber is an ex- 
tensive rigid body moving at changing angular frequency 
through an inhomogeneous sea of particles, which en- 
counter the perturbation in a periodic fashion. Yet, this 
series of experiments has demonstrated that the deceler- 
ation of the bar caused by dynamical friction from non- 
interacting halo particles scales with bar mass and halo 
density as predicted by eq. ||2Jl, which so closely resembles 
Chandrasekhar's formula (1). The naive scaling with ve- 
locity dispersion holds approximately for a general halo, 
and rather better for a self-similar halo, such as the SIS. 
The scaling with these parameters is an inevitable con- 
sequence of dynamical friction being second order effect 
caused by the interaction between the perturber and its 
own wake. 

6. EXCHANGES AT RESONANCES 

In this section, I present a more detailed analysis of one 
experiment in order to shed more light on resonant interac- 
tions, and the time-dependence of the frictional accelera- 
tion. The variation of Q{flb) for this "fiducial" experiment 
is shown by the dotted line that starts with f2t, — 0.5 in 
Fig.|2Ia) - i.e. the experiment with the lowest initial Qb, 
and the only one shown in this figure to start with coro- 
tation outside the bar. 

6.1. Adiabatic invariants 

The interaction between an orbit and a perturbing bar 
potential was discussed by Lynden-Bell (1979) in the anal- 
ogous case of disks. In that geometry, the unperturbed 
orbit of a particle in resonance closes in the frame that ro- 
tates with the perturbation, and Lynden-Bell pointed out 
that a star close to a resonance could be regarded as pur- 
suing a closed figure that precesses slowly relative to the 
major axis of the perturbation. 

When m — n, for the disk-like resonances, the unper- 
turbed orbit precesses at the rate f2m = (jifl^ + kil,r)/rn. 
In general, Qm ^ f^b, and the difference, — = is 
a slow frequency close to the resonance; i.e. \^sl^^\ ^ 1. 
Since the time taken for the star to complete an orbit 
round the closed figure is short (one or two radial peri- 
ods), the action J f = — kJ^/m, associated with this 
fast motion will be adiabatically invariant. On the other 
hand, the slow rate of precession of the orbit relative to 
the bar allows a large change to the "slow" action. 

The situation is more complicated for the DRR, since 
the orbit does not close in the bar frame and the adiabatic 
invariant is the total angular momentum, L (WK05). Con- 
servation of L might seemingly preclude any angular mo- 
mentum exchanges with the bar. However, this resonance 
is most important for highly inclined orbits, and while L 
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is conserved, the plane of the orbit, and therefore L^, is 
changed at the resonance; thus interactions at the DRR 
are able to slow the bar. 

Because even quite mild potential perturbations can 
produce non-linear orbital responses at resonances in this 
manner, TW84 examined the resonant terms more care- 
fully. They showed that second-order perturbation theory 
remains valid provided that particles make a "fast" pas- 
sage through the resonance. By this, they mean that the 
pattern speed of the bar is changing sufficiently rapidly 
that "large non-linear resonant perturbations do not have 
time to develop before the star has crossed the resonance." 
In effect, the star crosses the resonance in less than one 
libration period of the figure. Since the precession rate is 
slow, TW84 and Weinberg (1985) argued that fast pas- 
sage through the resonance is the correct assumption, and 
second-order perturbation theory should be valid. 

6.2. Individual orbits 

Figure IHl presents a direct test of this assumption in the 
fiducial run. The time evolution of 17;, is shown in panel 
(a), while panel (b) shows the evolution of the instanta- 
neous Lz for ten particles orbiting in the mid-plane. The 
orbits shown are those of test particles run in a reconstruc- 
tion of the time evolution of the total potential, which re- 
quires knowledge only of the bar phase at every step in the 
original simulation. All have the initial energy E = —0.38 
and angular momentum, ~ 0.4L,nax(£') but are spaced 
equally in initial phases. The short-period oscillations in 
Lz of each particle are forced as each particle moves in 
the bar field. Larger, and lasting, changes occur as parti- 
cles pass through a resonance. Individual orbits can either 
gain or lose, depending on their phase relative to the bar. 

The group of particles selected is typical of the many 
hundreds I have examined. Large changes in occur 
over the time interval t ~ 50 to f ~ 200 during which time 
Vtb decreases from ~ 0.45 to ~ 0.4. About half the parti- 
cles gain in angular momentum by a up to ~ 20% as they 
are caught in a corotation resonance, while the particles 
at other orbital phases lose similar amounts. These same 
particles also pass through an inner Lindblad resonance 
around t ^ 1100 where a majority gain Lz- Notice that 
friction is moderate when these particles pass through the 
CR, but that the late second blip in friction (Fig. oc- 
curs as the particles pass through the ILR. Many other 
in-plane orbits that I have examined also show two sepa- 
rate periods of angular momentum exchange with the bar. 

It should be noted that the more tightly bound the or- 
bit, the higher its intrinsic frequencies. Some of the most 
tightly bound orbits indeed pass through the ILR at early 
times, as stressed by WK05. However, these orbits are 
confined close to the center, where they couple weakly to 
the inner part of the perturbation, whose amplitude de- 
cays quadratically towards the center. Their consequent 
modest gains in angular momentum contribute little to the 
overall torque, which is dominated by the CR and OLR 
exchanges by larger orbits, as I show in ^6.31 The bar has 
already become uninterestingly slow by the time the larger 
orbits cause strong friction at the ILR. 

TW84 note that the frictional torque will be reduced, 
and the sign of the exchange between the particle and the 
perturbation difficult to predict, if the halo particles pass 
slowly through the resonance (as defined above). TW84 
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Fig. 8. — (a) The evolution of the pattern speed in the fiducial 
run. (b) The instantaneous angular momenta of ten particles as 
a function of time. The particles move in the mid-plane and all 
start with the same energy and angular momentum but have equally 
spaced phases. The quasi-periodic oscillations on each line are the 
forced response to the bar potential, the large changes, which are 
not modulated for the most part, are caused by passages through 
resonances, (c) An approximation to the radial action for the same 
orbits shown in (b). (d) The instantaneous z-component of angular 
momentum of a different ten particles that start on exactly polar 
orbits, but have the same energy and angular momentum as in (b) 
have equally spaced phases, (e) The total angular momentum of the 
particles shown in (d). 
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and Weinberg (1985) expected fast passage through the 
resonance, but offered no proof. Diagrams such as Fig- 
ure |Htb) aUow us to examine this key issue. The changes 
for each particle at CR are generaUy non-osciUatory, indi- 
cating that the particle is not librating in the resonance, 
and therefore the passage is fast. The story at the ILR is 
less clear, however, as the instantaneous of each par- 
ticle continues to oscillate right to the end. It is likely 
that these oscillations are simply the forced responses to 
the bar potential; their period seems to be about 100 time 
units, which is reasonable for precession of closed ellipti- 
cal figures with respect to the slowly-rotating bar. These 
oscillations aside, the larger changes in mostly, with 
perhaps a couple of exceptions, seem to occur within a 
single oscillation period. 

Thus Figure IHl^b), and hundreds of other orbits that I 
have examined, show that most passages through reso- 
nance can indeed be treated by the LBK approach. Since 
I have not examined every corner of phase space, this does 
not, of course, amount to a proof that the fast passage 
assumption is adequate, but it does suggest it is good for 
most orbits. 

The adiabatic invariance of the fast action is also borne 
out, as shown in Fig. |HIc). Determination of the ex- 
act radial and azimuthal actions in a rotating and slow- 
ing, non-axisymmetric potential would be very difficult, 
but approximate values are good enough for our purpose 
here. The instantaneous value of Lz — J^, and we adopt 
Jr ~ 0.5a^7r/T, where a is the semi-radial excursion of 
the particle and t, the radial half-period, is determined by 
the time between the most recent passages through peri- 
and apo-center for the particle. (Properly-defined actions 
would not exhibit short-period oscillations, and the non- 
smooth variations of Jr at late times reflect a change in 
the estimated value each time the radial direction of the 
particle reverses.) 

At the CR, fc = and therefore J/ = J^; it can be seen 
that the approximate radial action is almost unchanged 
during the changes to around t ^ 100. (Note the dif- 
ference in scales between panels b and c.) The opposite 
is true for the changes around t — 1100, where it is clear 
that particles that gain also lose Jr and vice-versa, as 
predicted by conservation of Jf when fc = — 1 for the ILR. 
In fact, AJr — —^l^zl'^ in these cases, confirming that 
resonance responsible for these large changes is indeed the 
ILR. 

Fig. ISfd) and (e) present a similar study of orbits that 
start with the same total angular momentum and en- 
ergy as in (b) and (c), but which are initially oriented 
perpendicular to the bar plane so that Lz = for all 
initially. These orbits pass through the DRR between 
300 ^ i <^ 400, and it can be seen that some gain 
and others lose. However, panel (e) shows that the total 
L of these orbits does not change during this period, con- 
firming that this resonance merely re-orients the plane of 
each orbit. (These same particles suffer changes to L at 
later times when they pass through other resonances.) 

6.3. Integrated changes 

It is interesting also to examine the changes to the dis- 
tribution of particles about the main resonances as fric- 
tion proceeds. Extracting this information from simula- 
tions is not completely straightforward because the an- 



gular momentum of a resonant particle on a nearly cir- 
cular orbit differs from that of another particle at the 
same resonance that has a highly eccentric orbit. HoUey- 
Bockelmann, et al. (2003) and WK05 draw contours of the 
net changes of angular momentum over a short time inter- 
val in (S, L/Lmax) space, where Lmax(-E') is the specific 
angular momentum of a circular orbit of specific energy 
E. In order to suppress shot noise, these figures require 
considerable smoothing, even when very large numbers of 
particles are employed. 

Perturbation theory (TW84, Weinberg 2004) indicates 
that friction depends on gradients of the DF w.r.t. both 
Jr and J^, although the coefficient of the radial action 
derivative is fc, which is zero at CR. Since Jf = Jr — kJ^/m, 
is adiabatically invariant, changes in Jr must be related to 
changes to J^, and we need examine the distribution with 
respect to a single action only. The natural choice is Lz, 
the quantity of greatest interest for friction. 

At any given resonance, the particles that undergo the 
largest changes have ilm — ^b- We define f{Jr, J 4,) to be 
the density of particles in the space of these two actions,^ 
and determine the mean value of / at fixed Om- The mean 
value is defined by the path integral 



f{Jr,J<p)dl 

Icdl - 



(10) 



with the path C being the locus of constant flm through 
(Jr, J^) space from = to = 0. In practice, the nu- 
merator is a kernel estimate from the values of flm com- 
puted from all the particles at each separate resonance. 
Since f2m is a single-valued function of the angular mo- 
mentum. Lies, of a circular orbit, I use this as the abscissa 



F{Lrcs) — 



J2t w(Lrcs - Lyc 



(11) 



where L^cs,i is the angular momentum of a circular orbit 
that has the same D,m as the j-th particle, w is a kernel 
function, and the denominator is unchanged. 

The average defined in eq. (|10|) is not a canonical trans- 
formation, and therefore does not preserve phase space 
volume. Thus we cannot expect the gradient of F(Lres) to 
be a completely reliable surrogate for the gradients in the 
DF - a Jacobian factor may alter the slope. But since F 
is defined as the average of / at fixed frequency difference 
from the resonance, we should expect its slope to give a 
general indication. 

Figure O shows the function F(Lrcs) at the ILR, CR, 
OLR and DRR at two separate times during the evolution 
of the fiducial model. A significant fraction of the particles 
in the simulation, of order 10%, contributes to the curve 
in each panel implying a rather low-level of shot noise. 

The strongest feature is for CR att — 200 (middle panel 
of top row), where F shows a clear peak just to the low-L 
side of the resonance, which is marked by the vertical line. 
Near to the resonance, particles cross corotation in both 
directions on horse-shoe orbits (Binney & Tremaine 1987, 
§7.5) as they gain and lose Lz- The negative gradient in 
F implies there are more particles on the low-L side, and 

^This is not the distribution function as usually defined, but is 
the DF integrated over all orbit phases and the inclination angle of 
the orbit plane 6. 
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Fig. 9. — The function _F(Lres) as defined in the text at the most important resonances at two difli'erent times in the fiducial simulation; 
the top row shows curves at t = 200, while the bottom row is for t = 1000. The vertical line shows location of the resonance at the moment 
illustrated, and the range of abscissae is set to show a region around the resonance containing ^ 10% of the halo particles. The vertical scale 
is linear from zero, but is otherwise arbitrary. 



therefore an excess of gainers over losers. This imbalance 
causes a net gain of angular momentum by the particles, 
and therefore a net loss by the bar - producing the required 
friction. 

If the pattern speed were to stay constant, the imbal- 
ance would tend to flatten the slope of F, and the dis- 
tribution of particles about the resonance would approach 
kinetic equilibrium in which there would be more nearly 
equal numbers of particles crossing in both directions. But 
as flp declines, the resonance keeps moving to larger Lrcs 
(frequency is decreasing function of angular momentum), 
and equilibrium is never established; instead, the density 
of particles about the dominant resonance (s) responsible 
for friction maintains a shoulder, or excess of particles, on 
the low-L side of the resonance, as shown in Fig. O 

The pronounced feature at CR at t = 200 is character- 
istic of the distribution over most of the evolution, and 
changes at CR dominate the torque. A similar, but less 
pronounced, feature can be seen at the OLR at t — 200, 
but its importance fades as the bar slows. The DRR also 
shows a definite, though weaker, feature at t = 200; note 
that the range of Lros seemingly overlaps that shown in the 
panel labeled CR, yet F{L^cs) is almost featureless near 
this value because the function is computed for a different 
resonance. The DRR feature persists for about the same 
duration as the CR, but is again too weak to be visible at 
the last time shown. 

The dominance of CR also fades in the latest stages 
when the bar speed has more than halved. A very mild 
flattening of F at the ILR is possibly present at t = 1000, 
and is the only indication I could find in plots of F{Ly-cs) of 
changes to i^(ircs) at this resonance. However, the weak- 
ness of features in F(Li.os) at this resonance does not imply 
that no changes occur; indeed Fig. |S1 shows that substan- 
tial interactions are taking place at the ILR around this 
time. The weakness of this feature at the ILR may indicate 
that exchanges at this resonance are not self-reinforcing, 
in contrast to those at other major resonances. Alterna- 
tively, features may be suppressed by two effects: the most 



important is that the instantaneous angular momenta of 
particles that have passed through the ILR oscillate by 
amplitudes of typically 0.05 (see Fig. IS];), which may be 
sufficient to wash out intrinsic features. Second, conser- 
vation of Jf at this resonance implies that changes in 
must be reflected in changes of the opposite sign in L2, 
leading to partial cancellation in the net change in L^. 

7. TIME DEPENDENCE 

Weinberg (2004) finds that the LBK torque formula 
from TW84 does not predict the correct time evolution of 
the bar pattern speed; he was able to obtain approximate 
agreement with his experimental results only by taking the 
history of the bar perturbation into account. His formula 
for the torque in the time-dependent regime is also of the 
form eq. JSjl, but with a function Q that now depends on 
the entire history of the perturbation. 

Since a theoretical prediction invites comparison, I have 
tried experiments similar to that reported by Weinberg. 
The solid curve in Figure [TUI shows the pattern speed evo- 
lution using my code with the bar and NFW halo (see 
Appendix) that Weinberg employed for his first test. The 
axes have been scaled to the units adopted by Weinberg 
and I have reproduced the curves from his figure (see also 
Fig. 16 of WK05); his time-dependent, linear theory pre- 
diction is shown by the dashed curve and Weinberg's own 
simulation is shown by the dot-dashed curve. Both curves 
differ from each other and from my result, although the 
differences are minor. 

The small differences from both Weinberg's theoretical 
prediction and from his simulation do appear to be sig- 
nificant, however. I have checked that my result does not 
depend on any numerical parameters; increasing the num- 
ber of particles to lOM, halving the time step, increasing 
the outer truncation radius by 50%, and other tests all 
yielded results that differ from the solid curve shown by 
scarcely more than the thickness of the line. I have also 
verified that total angular momentum is conserved to 0.3% 
of the small total possessed by the bar at the start. 
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Fig. 10. — The dashed hne reproduces Weinberg's (2004) predic- 
tion from perturbation theory for the pattern speed evolution, while 
the dot-dash line shows his Af-body simulation result; the model he 
used is described in the Appendix. The solid line shows the evolu- 
tion in my experiment for the same case. See j2|for details. 

Weinberg's prediction requires sophisticated numerical 
treatment of a stiff system of coupled equations, and dif- 
ferences from my simulation could possibly result from nu- 
merical inadequacies in his prediction. It is harder to un- 
derstand why Weinberg's simulation differs from mine, as 
I have found my result to be so robust. 

Returning to Fig. [2Ia), the various dotted curves indi- 
cate a strong, and long-lasting, dependence on the initial 
pattern speed. I have examined orbits and F for the dra- 
matically stronger friction exhibited by the curve shown 
by the solid line in this figure, in the same manner as for 
the fiducial model reported in I find that the stronger 
drag at a given f2h when the pattern speed is started from a 
higher value arises from a steeper local gradient in F{Ly-cs) 
because the shoulder is more pronounced than that shown 
in Fig. |51 The resonance lies in similar position relative 
to the peak, but the higher peak has steeper slopes. The 
shoulder is more pronounced because the higher starting 
speed enables more particles at higher frequencies to be 
caught in the CR, and to experience much larger changes 
in Lz (some gain by up to factors of 5). 

7.1. Discussion of time- dependence 

In hindsight, the importance of time-dependence is en- 
tirely reasonable. It is clear from the distribution of par- 
ticles about the CR and OLR shown in Fig. El that the 
local slope of F depends on the past history of resonant 
interactions. As it develops, friction sculptures the DF in 
the vicinity of the dominant resonance, which in turn af- 
fects the strength of the subsequent frictional drag. Thus 
the limit of a steady bar rotating at fixed fif,, which was 
adopted by TW84, does not capture the true story; Wein- 
berg's (2004) re-derivation using Laplace transforms is es- 
sential. 

The very strong dependence of the acceleration on the 
starting angular frequency reported in Fig.|2a) is reflected 



in the magnitude of the shoulder, and the local gradient 
in F at the resonance. The curves in that figure seem to 
converge for smaller values of fit,, suggesting somewhat less 
extreme dependence on past history in the physically more 
reasonable regime where corotation is beyond the end of 
the bar. 

Figs, in -[3 show that the scaling of the frictional force 
with the main parameters of eq. (j3Jl still holds, at least ap- 
proximately, despite the complication of time-dependence, 
as Weinberg also notes. Note that the perturbation has the 
same initial Oh and turn-on rule for the bar mass and halo 
density tests, so that the time-dependence of the function 
Q did not change. However, the initial bar pattern speed 
had to be scaled in the experiments with different veloc- 
ity dispersion, making the reasonable agreement shown in 
Fig. EKb) all the more reassuring. 

8. CONVERGENCE TEST 

As noted earlier, halo particles gain angular momen- 
tum, on average, at resonances. Even though each may 
pass through the resonance rapidly, we will not obtain 
the smooth torque expected in the continuum limit un- 
less there are many particles, densely spread over a broad 
range of frequencies. The required number depends on the 
width of the resonance in frequency space: If the bar had 
a fixed or very slowly changing pattern speed, the reso- 
nances would be sharp and the simulations would indeed 
require a very large TV (Weinberg & Katz 2002; HoUey- 
Bockelmann, et al. 2003; WK05), but the forcing frequency 
is broadened by its decreasing angular speed and we can 
hope that the discreteness of the particle distribution will 
become insignificant for some attainable N. 

Figure ITTI shows that the pattern speed evolution in the 
Hernquist halo appears to converge to a smooth function 
as N increases. Merely 10^ particles are sufficient to obtain 
qualitatively correct behavior for the 5% mass bar, but al- 
most one hundred times larger TV is needed for comparable 
numerical quality with a bar of one tenth this mass. Phase 
space needs to be populated more densely as the mass of 
the bar decreases, in order that a large enough number of 
particles are passing through the resonance at any time to 
yield a smoothly varying force. Both the effective potential 
of the resonance is weaker, and the Lorentzian frequency 
width of the resonance is reduced by the slower braking 
rate. As both factors vary with the mass of the perturber, 
it seems reasonable that the number of particles needed 
to obtain smoothly varying friction should increase faster 
than . However, one million is adequate even for a 
low-mass bar. 

The smoothness of the braking rate in Figs.|21&EI and 
for large N in Fig. 1111 is evidence that the resonant ex- 
changes which give rise to the force are dense enough at 
most frequencies to produce a smooth torque. The fric- 
tional force may become erratic at very low pattern speed 
(Fig. [2Ji) because phase space is not populated densely 
enough as \ flb\ decreases, but this is a physically uninter- 
esting regime. 

These empirical results conflict with the recent claims by 
WK05 that a much larger number of particles is needed to 
obtain the correct result. However, their analysis does not 
take proper account of the changing pattern speed, which 
is so important. The time dependence of the frictional 
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Fig. 11. — The bar pattern speed as a function of time for exper- 
iments with differing mass bars, each for four separate values of N. 
The indicated bar masses are fractions of Mj, and the length of the 
time axis is inversely proportional to the bar mass in each panel. 
The line styles indicate the number of particles: A'^ = 10* - dotted, 
AT = 10^ - dashed, AT = lO'' - dot-dashed, and = 10^ - soHd. 

force arises from preceding changes to the DF, and the 
shoulders in Fig. El indicate directly just how broad the 
resonance is. The feature can develop and behave as it 
does only if particles with angular momenta as far away 
as the width of the shoulder from the precise resonance are 
being affected by the resonance. Since ~ 10% of all halo 
particles contribute to each panel, it appears that ^1% 
of all halo particles take part in exchanges at a dominant 
resonance at any moment. The convergence tests in Fig. 1111 
confirm that stochastic behavior when resonances are too 
sparsely populated (what WK05 describe as "coverage" ) is 
brought under control at comparatively low-iV, when the 
pattern speed changes at realistic rates. 

It should be noted that these experiments reveal only 
the number of particles needed to populate phase space 
densely enough for the frictional force to approximate the 
continuum limit in the absence of self-gravity. The issue 
of coUisional relaxation also needs to be addressed when 
the particles interact with each other, which may require 
more particles, as I discuss in ^19.31 

9. EFFECT OF SELF-GRAVITY IN THE HALO 

All experiments reported so far treat the halo as a non- 
interacting population of test particles, which move in a 
fixed potential well and respond only to the perturber, 
as in the usual treatment of dynamical friction {e.g. BT). 
Collective effects in a collisionless halo could also be im- 
portant, however, because the wake itself contributes to 
the non-axisymmetric density affecting the orbits of the 
halo particles. Hernquist & Weinberg (1992), Weinberg & 
Katz (2002), and Sellwood (2003) have already reported 




Fig. 12. — The acceleration of the bar when self-gravity of the halo 
response is included. The solid curve shows the behavior when all 
terms < ( < 4 are included, the dashed curve shows the situation 
when only the 1 = 2 terms are employed. The dotted curve shows 
the corresponding case with no self-gravity, the fiducial model, re- 
produced from Fig. l2la'>. 

some experiments with rigid bars that take this into ac- 
count, but the emphasis in those papers was on the change 
to the halo density profile. 

Here I study the effect of self-gravity on the frictional 
force, while still employing an imposed bar. I break the 
discussion of collective effects into two levels of compli- 
cation: when self-gravity includes the monopole terms, 
the radial mass profile of the halo could possibly change 
substantially over time, as claimed by Weinberg & Katz 
(2002). I therefore begin by including only the quadrupole 
field of the response density, before describing how other 
terms affect the evolution. Note that all the experi- 
ments described in this paper are perturbed with the non- 
axisymmctric field of a rigid bar. Fully self-consistent sim- 
ulations, with the bar also made of responsive particles, are 
reported in other work {e.g. Papers II & III). 

I compute the self-gravity of the halo using the method 
described by McGlynn (1984) and Sellwood (2003, Ap- 
pendix A). Briefly I use a 1-D spherical grid, with an ex- 
pansion of the gravitational field up to order /max in surface 
harmonics on each radial shell of the spherical grid. 

9.1. Quadrupole field of the response 

The dashed curve in Figure IT^ shows the bar accelera- 
tion when only the quadrupole field of the halo response 
density is included; the dotted curve, reproduced from 
Fig. Ola), shows the behavior with no self-gravity. It is 
clear that the self-gravity term increases the drag slightly 
when the pattern speed is slow, but greatly diminishes it 
when the bar has an artificially high pattern speed. 

This behavior is a consequence of the phase of the halo 
response. Including the self-gravity of the response weak- 
ens the net torque on the bar when the response is more 
than 45° out of phase with the imposed bar; the phase 
lag when self-gravity is included is somewhat similar to 
that shown in Fig. Elb) , but the response remains more 
nearly perpendicular for longer. As the bar slows, the re- 
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Fig. 13. — The evolution of the pattern speed of a 2% mass bar 
which starts with corotation at the end of the bar (a = rg), for 
experiments with different A^. Unlike for Fig. 1111 these experiments 
include all low-order self-gravity terms of the halo particles (Z < 4) . 

sponse gradually shifts towards alignment with the bar, 
thereby augmenting the quadrupole field of the perturba- 
tion and causing modestly increased friction. Note that 
friction peaks at fib — 0.5 when corotation is at the end of 
the bar, implying that self-gravity always enhances friction 
in the physically relevant regime. 

9.2. Including the monopole and other terms 

The solid line in Fig. II 21 shows the effect of adding more 
terms to the self-gravity of the halo response; in this case 
all / < 4 terms. The addition of these extra self-gravity 
terms increases friction somewhat at most angular speeds 
over that obtained when only I = 2 terms are used. 

9.3. Convergence test 

FigureElshows the effect of changing the particle num- 
ber, N . I use the Hernquist halo, a bar with Mi, — 0.02M/i, 
and a more realistic initial 17^. This test includes alH < 4 
terms of the self-gravity of the halo density response and 
thus some coUisional relaxation must be present. 

As without self-gravity, the time variation of the pattern 
speed becomes smoother as TV increases. The evolution 
is closely similar for the two cases with N > 10^, but 
friction is clearly overestimated when N = and slightly 
so when N — 10^. The minimal differences between the 
results for the experiments with N = 10^ and N — 10^ 
indicate that IM particles is a sufficient number to capture 
the correct physics for this case with a low-mass bar. 

The 2% mass bar used in Fig. ^]is the same as for the 
second most strongly braked case in Fig. ^] Comparison 
of these two convergence tests indicates that the low-iV 
models depart more strongly from the high-iV results when 
self-gravity is included. However, the number of particles 
needed for convergence is not increased dramatically by 
self-gravity. Thus, Weinberg & Katz (2002) and HoUey- 
Bockelmann, et al. (2003) are correct that the reduction in 
orbit quality caused by numerical noise in self-gravitating 



halos does affect the friction force, but the torque con- 
verges for reasonably accessible numbers of particles; IM 
appears to be plenty in this problem, as Sellwood (2003) 
concluded from fully self-consistent experiments. 

Note also the complete absence of numerical noise in 
the higher- iV experiments shown in Figs. & El The 
convergence is truly impressive; curves for which N dif- 
fers by factors of 10 overlay almost perfectly, and there is 
little evidence of differences due to shot noise. ^ It seems 
inconceivable that such impressive convergence could be 
achieved if Weinberg's arguments {e.g. Weinberg & Katz 
2002 and subsequent papers) for the much higher number 
of particles needed were correct. 

In this vein, it should also be noted that my earlier 
work to reproduce the normal modes of disks {e.g. Sell- 
wood 1983; Earn & Sellwood 1995) already demonstrated 
that simulations with modest numbers of particles could 
capture global instabilities, which are driven by resonant 
interactions between particles and rotating potential per- 
turbations. The eigenfrequency predicted from linear the- 
ory can be reproduced to a precision of a percent or two 
in simulations with a few tens of thousands of particles. 
Admittedly these results were for a dynamical instability 
in a thin disk, whereas halo friction is a secular effect in 
3-D. However razor-thin disks with random motion are 
scarcely more complicated than a spherical model, where 
each orbit is confined to its own plane. Furthermore, since 
the evolution in both cases is determined by the resonant 
terms, they have similar requirements for phase coverage 
of particles at the resonance, the ability of the resonance 
to persist for a large number of dynamical times, etc., and 
there is no clear reason why particle number need be or- 
ders of magnitude higher for secular evolution than for 
instability. 

9.4. Claimed counter-example 

HoUey-Bockelmann, et al. (2003) report fully self-con- 
sistent experiments for which the behavior with N = IM 
differs from similar experiments with larger TV. However, 
the difference does not necessarily support their conclusion 
that IM particles is inadequate for bar-halo interaction. 
They show (their Fig. 20) that the quadrupole field in the 
IM particle experiment is weaker, so the halo torque must 
also be weaker, which implies that the lower halo particle 
number could still be adequate. Both Weinberg & Katz 
(2002) and Sellwood (2003) report that smaller N leads 
to more rapid angular momentum transfer, as also shown 
in Fig. 1131 making the inadequate- iV interpretation in this 
case still less plausible. HoUey-Bockelmann, et al. do in- 
deed have a result that is iV-dependent, but the effect of 
lower particle number is to decrease the strength of the 
m = 2 distortion in the disk, and is therefore a conse- 
quence of disk dynamics and not halo friction. It is likely 
that a higher noise level in their lower- TV disk interferes 
with their bar triggering mechanism. Specifically, they 
start with a bar-unstable disk and also apply a transient 
tidal field to trigger the bar; the smaller the number of 
particles, the larger the initial amplitude of the intrinsic 

•^This aspect is, in part, due to the careful selection of particles 
from the DF for the initial set-up (see Debattista & Sellwood 2000 
Appendix B); random selection of particles inevitably leaves \/N- 
type variations in the density of particles as a function of the integrals 
and does not, therefore, achieve quite such impressive convergence. 
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instability, which will generally interfere with the applied 
tidal field leading to a weaker bar, unless the phases of 
the two bar-forming mechanisms happened to be nearly 
aligned. 

10. CONCLUSIONS 

I have shown that dynamical friction between a rotat- 
ing, imposed bar and a halo scales in a very similar man- 
ner to that predicted by Chandrasekhar's formula. The 
experiments reported in §21 & El are highly idealized and 
employed the simplest possible system that could capture 
dynamical friction on a bar rotating in a halo. The re- 
tarding acceleration depends on the angular speed of the 
bar, scales linearly with its mass {i.e. the strength of its 
quadrupole field) and with the background density. It 
also scales roughly inversely as the square of the velocity 
dispersion of the background through which it moves, al- 
though this scaling is more approximate unless the model 
is self-similar. This result is the analog for a bar of that 
obtained by Lin & Tremaine (1983) from a similar study 
with orbiting satellites. 

Even though the physical situation is quite different 
from rectilinear motion through an infinite, uniform back- 
ground, and complicated by the existence of resonances 
and bar turn-on issues, the Chandrasekhar scaling still 
holds. It holds, because dynamical friction is fundamen- 
tally a second order effect that arises from the interaction 
of a perturber with its own wake. The strength of the in- 
teraction does depend on the details, but the parameter 
scaling cannot. 

In all three halo types employed here, friction is weak 
when the bar pattern speed is so high that corotation is 
well inside the bar. As the bar slows, the frictional drag 
grows at first but generally peaks before corotation reaches 
the bar end. Thus, in the physically interesting regime, 
where corotation is beyond the end of the bar, friction 
always decreases as the bar slows in these non-rotating 
halos. 

Friction is dominated by the quadrupole field of the 
bar both because the quadrupole is the dominant poten- 
tial component of the bar, but also because the higher- 
order resonances that are associated with the higher ex- 
pansion terms couple less strongly to the particles. Since 
the lowest-order, non-axisymmetric component of the bar 
field dominates, friction can be captured adequately in 
simulations with even quite low spatial resolution, pro- 
vided enough particles are employed. 

Tremaine & Weinberg (1984) first demonstrated that 
friction arises as the perturbation sweeps across resonances 
with the particles, a point stressed in recent work by Wein- 
berg and his co-workers and by Athanassoula (2003). I 
have shown that the pattern speed of the bar changes 
sufficiently rapidly that a halo particle generally passes 
through the resonance without any complicated non-linear 
trapping, as expected by TW84. The forced responses of 
the halo particles change from ant i- alignment to alignment 
with the bar as the pattern speed crosses a resonance, 
which happens smoothly because the resonance is broad- 
ened by the changing pattern speed. The net effect is to 
produce a global density response that lags the bar, caus- 
ing the frictional drag. The lag angle between the bar and 
the halo response varies with the friction force, as shown 
in Fig.|21[b); it is close to 45° when friction is strong, but 



the halo response is closely aligned with the bar when the 
bar is slow, and almost orthogonal when the bar is unrea- 
sonably fast. 

Since one of the two actions associated with the un- 
perturbed motion is adiabatically invariant, exchanges at 
resonances can be examined as a function of a single an- 
gular momentum-like variable, which I denote Lrcs- The 
density of particles is generally a decreasing function Lres, 
although previous angular momentum exchanges with the 
bar cause a local shoulder to develop at the most im- 
portant resonances. The distribution with L^es evolves 
most strongly at corotation at physically interesting pat- 
tern speeds, but smaller changes are detectable at other 
resonances before the bar has slowed much. 

The sculpturing of the particle distribution about the 
principal resonances by previous friction seems to be the 
reason for the strong dependence of friction on the previous 
evolution, discovered by Weinberg (2004). The magnitude 
of the friction force is determined by the local gradient in 
F, which differs from that in the initial model because pre- 
vious exchanges with the perturbation have rearranged the 
particle distribution. A possibly similar effect may hap- 
pen in proto-planetary disks ( Artymowicz 2004) , although 
self-reinforcing responses in that context have been found 
only at corotation, whereas I have found them at other 
resonances also. 

In the absence of self-gravity, the number of particles 
needed to obtain a smoothly varying frictional force is 
quite modest, unless the bar is very weak. The main re- 
quirement here is that the broadened resonances caused 
by the time-varying pattern speed should overlap many 
particles in order to obtain the correct density response. 
The pattern speed changes more slowly for weaker bars, 
decreasing the frequency width of the perturbation and 
consequently raising the particle number needed to obtain 
a smoothly varying force. The number of particles needed 
in this regime appears to rise more steeply than the inverse 
of the quadrupole field strength. 

Weinberg & Katz (2002) and WK05 argue that coUi- 
sional relaxation in simulations with self-gravity should 
increase the number of particles required to approach 
the continuum limit, because potential fiuctuations aris- 
ing from Poisson noise affect the orbital behavior. The 
lowest-energy orbits, which they find are most delicate, 
should make a negligible contribution to the torque. I have 
found a small reduction in force quality when self-gravity is 
included, consistent with the prediction by these authors, 
but quite modest particle numbers are needed to bring this 
problem under control. Thus the frictional torque can be 
simulated accurately with standard algorithms and readily 
accessible computers. 

All experiments reported here employ an imposed bar 
potential in order to examine the dependence of the halo 
response on the bar parameters. This simplifying approx- 
imation has many obvious disadvantages, as noted in ^ 
Fully self-consistent simulations, such as those to be re- 
ported in later papers in this series, are the only way to 
ensure that the bar forms with a dynamically realistic am- 
plitude and pattern speed and responds to friction in the 
appropriate way. 
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potential when r < rg are almost the same, while the 
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He chose a large, strong, heavy bar, with semi-major axis 
equal to = 1/15, axis ratios a : & : c = 1 : 0.2 : 0.05, and 
a mass equal to half the halo mass interior to this radius, 
or 0.0526 of the virial mass. The initial pattern speed of 
the bar is 18.84 in his units, which places corotation at the 
bar end. The bar amplitude varies as f{t) = [1 + erf (4t — 
2)]/2, where t is the time in Gyr for a halo scaled to the 
Milky Way; when ~ 20 kpc, Weinberg's time unit is 
- 1.4 Gyr. 
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APPENDIX 



The NFW (Navarro, Frenk & White 1997) halo density 
profile is 



/Onfw 



Part 



r{r + r^)^ 



(1) 



where pa is a scale density and is a scale length. The 
concentration parameter c is defined as the ratio of the 
"virial radius" to Vs', at the virial radius, the average en- 
closed density is 200pcrit = 200 x 3i?^/(87rG), with 
being Hubble's constant. Eddington's formula (Binney & 
Tremaine 1987, eq. 4-140b) yields an isotropic distribution 
function that is positive everywhere, but which requires 
some care to evaluate for the most bound energies. 

Note that the NFW profile has a more slowly declin- 
ing outer density gradient than does the Hernquist profile 
(eq. but is otherwise closely similar. The density and 



